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Abstract 

Measures of nonlinearity and complexity, and in particular the study of Lyapunov exponents, have been increasingly used to 
characterize dynamical properties of a wide range of biological nonlinear systems, including cardiovascular control. In this 
work, we present a novel methodology able to effectively estimate the Lyapunov spectrum of a series of stochastic events 
in an instantaneous fashion. The paradigm relies on a novel point-process high-order nonlinear model of the event series 
dynamics. The long-term information is taken into account by expanding the linear, quadratic, and cubic Wiener-Volterra 
kernels with the orthonormal Laguerre basis functions. Applications to synthetic data such as the Henon map and Rossler 
attractor, as well as two experimental heartbeat interval datasets (i.e., healthy subjects undergoing postural changes and 
patients with severe cardiac heart failure), focus on estimation and tracking of the Instantaneous Dominant Lyapunov 
Exponent (IDLE). The novel cardiovascular assessment demonstrates that our method is able to effectively and 
instantaneously track the nonlinear autonomic control dynamics, allowing for complexity variability estimations. 
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Introduction 

Hearth contractions are regarded by many scientists as the 
foremost example of a physiological system showing predomi- 
nantly nonlinear behavior, mainly generated through integration 
of multiple neural signaling at the level of the sinoatrial node [1]. 
Accordingly, the fluctuations in the interval between consecutive 
heartbeats have been widely investigated as output of a nonlinear 
system revealing and quantifying the complexity of cardiovascular 
control [2-27]. 

Among all nonlinearity and complexity measures, Lyapunov 
exponents (LEs) have been proven to provide an important 
mathematical tool in characterizing dynamical properties of a 
nonlinear system [28] . LEs were first defined by Lyapunov [29] in 
order to study the stability of non-stationary solutions of ordinary 
differential equations and for more than fifty years they have been 
extensively studied in many disciplines [30-33]. Specifically, they 
refer to the average exponential rates of divergence or conver- 
gence of neighboring trajectories in the system phase space. In 
fact, for a system whose characteristic equations are known, there 
is a straightforward technique for computing the whole Lyapunov 
spectrum [28]. Several methods for a reliable data-driven LEs 
estimation, even in short time data records, have been also 
proposed [34,35]. In a deterministic nonlinear system with no 
stochastic inputs, a positive EE reflects sensitive dependence to 
initial conditions and can be taken as a definition of a chaotic 
system [36] . Nevertheless, a small amount of noise in a limit cycle 



oscillation could yield a positive EE if the trajectory has regions 
with large slopes. In chaotic systems [37], stochasticity does not 
play a crucial role and their dynamics are highly dependent on the 
initial condition. Although it is straightforward to consider chaotic 
mathematical systems in which stochastic inputs are suppressed, in 
actual applications especially related to physiological systems, it is 
not possible to eliminate such inputs making the chaos assessment 
simply unreliable [3]. Stationary aperiodic behavior, in fact, can 
also arise in linear or nonlinear stochastic systems. In light of these 
considerations, as this work deals with (instantaneous) LEs 
estimation with applications on heartbeat dynamics, we do not 
address the issue related to the chaotic behavior of heart rate 
variability (HRV). 

Relying on the approach suggested by Chon et al. [38] and, 
later, by Armoundas et al. [39], we consider the cardiovascular 
system both chaotic and stochastic. This concept is in agreement 
with current physiological knowledge, since healthy HRV 
dynamics can be considered/ modeled as the output of a nonlinear 
deterministic system (the pacemaker cells of sinus node) being 
forced by a high-dimensional input (the activity in the nerves 
innervating the sinus node). Accordingly, we model the heartbeat 
nonlinear dynamics as a third-order Nonlinear Autoregressive 
(NAR) model embedded in a point process probabilistic frame- 
work. Such statistical approach allows us to estimate the LEs in an 
instantaneous fashion by fitting the model to the observed data 
and applying the Fast Orthogonal Search (FOS) algorithm [40]. 
Point-process theory has been widely recognized as an excellent 
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mathematical tool to characterize the probabilistic generative 
mechanism of the heartbeat at each moment in time [41]. In the 
considered model, the intrinsically discrete, unevenly spaced 
heartbeat intervals are represented by a physiologically-plausible 
inverse-gaussian (IG) distribution. Defining the first and second- 
order moments of the IG distribution as function of the past 
heartbeat intervals (i.e., the RR intervals), it is possible to obtain an 
effective prediction of the next heartbeat event together with an 
accurate assessment of instantaneous indices of cardiovascular 
control. In previous studies [19,41,42], we demonstrated how to 
estimate heartbeat dynamics even in short recordings under 
nonlinear and non-stationary conditions using Wiener- Volterra 
theory for nonlinear systems identification. 

Introduction to the Instantaneous Dominant Lyapunov 
Exponent 

In this work, the IG mean is modeled as a third-order nonlinear 
function of the past RR intervals. In order to perform an effective 
parameter estimation and retain all the historical information of 
events, the cubic NAR kernels of the Wiener- Volterra series are 
expanded using the Laguerre bases [43] leading to the definition of 
cubic Nonlinear Autoregressive Laguerre (NARL) model. Of note, 
the NARL definition includes an infinite regression of the past 
events with a parsimonious use of model parameters. From the 
NARL point-process model and Fast Orthogonal Search algo- 
rithm, we are able to estimate the complete LE spectrum at each 
moment in time, providing novel information concerning the 
complexity dynamics and its variability. Of note, to the best of our 
knowledge, complexity variability measures have never been 
estimated from instantaneous indices of complexity and could 
open novel perspectives on the assessment of discrete stochastic 
physiological systems. We present two applications on synthetic 
datasets (the Henon map and Rossler attractor) and two 
experimental applications portraying the crucial role of the 
instantaneous Lyapunov Exponents in assessing autonomic 
changes in humans (ten heathy subjects undergoing postural 
changes, and fourteen patients with severe heart failure), focusing 
our attention on the Instantaneous Dominant Lyapunov Exponent 
(IDLE, A), which is the first exponent of the Lyapunov spectrum. 
Preliminary results of these analyses have been presented in [44— 
46]. We start with a detailed, exhaustive presentation of our 
methodological framework through the following "Materials and 
Methods" section. 
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where 7it = (uj,RRj,RRj-\,...,RRj-M+ \) is the history of the 
point process, i(t) is the vector of the time-varying parameters, 
IJ.jiji(t,7it,i(ty) represents the first-moment statistic (mean) of the 
distribution, and ^(){t) = 0>0 denotes the shape parameter of the 
IG distribution (as 6/ /i^go, the IG distribution becomes more like 
a Gaussian distribution). As f(t\Ht,^(t)) indicates the probability 
of having a beat at time t given that a previous beat has occurred 
at Uj, iiijiji{t,7it,i{ty) can be interpreted as signifying the average 
(or expected) waiting time before the next beat. We can also 
estimate the second-moment statistic (variance) of the IG 
distribution as o\j^{t) = fi\j^{t) / 0 . The use of an IG distribution 
to characterize the R-R intervals occurrences is motivated by the 
fact that if the rise of the membrane potential to a threshold 
initiating the cardiac contraction is modeled as a Wiener process 
with drift, then the probability density of the times between 
threshold crossings (the RR intervals) is indeed the inverse 
Gaussian distribution [41,49]. As a matter of fact, when compared 
with other distributions, the IG model achieves the overall best fit 
in terms of smaller KS distance [42] . The instantaneous RR mean, 
fijiji(t,7it,^(t)), can be modeled as a generic function of the past 
RR values fij^j^(t,ntAt))=g(^^j-U^^j-2,-,RRj-h), where 
RRj^}^ denotes the previous k^^ R— R interval occurred prior to 
the present time t. 

Nonlinear Modeling of History Dependence 

The general Nonlinear Autoregressive (NAR) formulation can 
be written as: 



y{n) = ¥(y(n - l),3;(n - 2),. . . ,y(n - M)) + e{n) 



(2) 



where e(n) are independent, identically distributed (i.i.d.) Gaussian 
random variables. The expected value of 3;(^) can be written as a 
Taylor expansion under the hypothesis of being infinitely 
difierentiable in a neighborhood of the working point: 



Materials and Methods 

Point-Process Models of Heartbeat Dynamics 

A random point process is a stochastic process whose elements 
are point patterns specified as a locally finite counting measure 
[47]. Considering the R-waves detected from the Electrocardio- 
gram (EGG) as such events, point process theory can be used to 
characterize their probability of occurrence [41,42,48]. Mathe- 
matically, in the time domain, a simple 1 -dimension point process 
consists of series of timestamps marking the occurrence times 
^ G [0,oo) of the random events. Given a set of R-wave events 
{uj}j^l, let RRj = Uj — Uj-i>0 denote the R-R interval, or 
equivalently, the waiting time until the next R-wave event. 
Assuming history dependence, the probability distribution of the 
waiting time t — Uj until the next R-wave event, where Uj denotes 
the previous R-wave event occurred before time /, follows an 
inverse Gaussian (IG) model: 



M 

E\y{n)\=yQ + ^yi(i)y{n-i) 

i=\ 

(3) 

00 M M K ^ ^ 

+ ^ ^ . . . ^ y^(/i , . . . ,iK)Tl y{n - ij) 

Due to the autoregressive structure of eq. 3, the system can be 
identified with only exact knowledge of the output data and with 
only few assumptions on the input data. In practice, this series will 
obviously be truncated at some small, finite value. Here, we 
represent the nonlinear cardiovascular system by taking into 
account up to the cubic nonlinear terms, i.e. yg? 7105 Yii^J)^ 
73 (/,/,/(:). Thus, the model can be written as: 
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M 

E\y{n)\=yQ + ^y\{i)y{n-i) + 

M M 

y2(ij)y(n - i)y{n -j) + (4) 

/=ly=l 
M M M 

XI XI XI y^^^J^^)y^^ -j)y{n - k) 

i=\ j=\ k=\ 

where the quadratic and the cubic terms, yiiQ) and y^{ij,k), are 
assumed to be permutation invariant. This choice of a third order 
NAR system further gives robustness against the presence of 
measurement noise in the data [38]. 

Laguerre Expansion of the terms. An important 
practical hmitation in modeling high-order nonlinearities using 
the model in eq. 4 is the high number of parameters needed to 
sufficiently fit the observed data. An advocated approach to solve 
such a limitation has been proposed by using the Laguerre 
functions [43,50,51]. Let us define the y^^-order discrete time 
orthonormal Laguerre function: 

0/«) = o<¥(l-a)I YS-\y(^^ (^Qoe'-'(l-a)', 

where a is the discrete- time Laguerre parameter (0 < a < 1) which 
determines the rate of exponential asymptotic decline of these 
functions, and n>0. Given the Laguerre function, ^j{n), and the 
signal, y{n), the /^-order Laguerre filter output is: 

GO 

lj{n)=Y, ^j{i)y{n-i-\) (5) 

The computation of the Laguerre Filter output can be simplified 
by using the following recursive relation [43] : 

h{n) = ^/oLk{n-\) + V\^y{n - 1) (6) 



lj{n) = ^lj{n-\) + ^ Ij- 1 {n) + (7) 

V^lj.xin-\)J>\ (8) 

Since the {0/(0} form a complete orthonormal set in functional 
space £2, we can write [52]: 

yo=^o (9) 



p 

yi(()=^gi(r)m (10) 



Q Q 

y20V) = XX S2{r,s)^XMs (J) if 
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yi(ij,k)= E E giir,s,tWi)mm) if i^J ■ (12) 

r=0 s=0 t=0 

Here go, gi(r), g3(r,s) and g3(r,s,t) are the Laguerre 
coefficients. Using eq. 5 and eqs. 9-12, the model in eq. 4 for 
the instantaneous RR mean becomes: 

p 
i = 0 

Q Q 

XX^2(^'^',owo/xo+ (13) 

/=0 7=0 

XXX^3(v,M/.(0/X0 4(0. 

/ = 0 7 = 0 k = 0 

hereinafter called Nonlinear Autoregressive with Laguerre expan- 
sion (NARL) model. Here, the Laguerre filter outputs are: 

m 

li(t)= ^Un)(^^N(t)-n-^^m)-n-l) (14) 

n = l 

with N(t)= lim^^;- N(t)= m3x{k : Uk<t} sls a. left continuous 
function. 

The coefficients go,{gi(i)}, {giiij)}, and {g3(ij,k)} corre- 
spond to the time-varying zero-, first-, second-, and third-order 
NARL coefficients, respectively. When a = 0 the filter output 
becomes lj(n) = {—iyy{n—j—l) and the NARL model corre- 
sponds, apart for the sign, to the finite NAR model in eq. 4, 
whereas for a 7^ 0 the instantaneous RR mean in terms of NAR 
equations is theoretically defined as follows: 

GO 

H^^(t,ntm) = ^^Nit)-i +>^o+ X ^'i^'"'^) ARR,- + 

00 00 

E E y^iiJ^t) ARR,- ARR,- + (15) 

;=i 7=1 

00 00 00 

+ X X X ^3(^,^,0 ARR. ARR.ARRyt 

i=i j=i k=i 

where ARR/^ = (RR - ~ ^^^^ N(t) - /? - 1 ) • The autoregressive 
model is expressed in terms of the derivative RR series, rather 
than the original RR series, in order to allow the model to cope 
with highly non-stationary regimes [53]. 

Moreover, as iij^j^{t,7it,^{t)) is defined in a continuous-time 
fashion, we can obtain an instantaneous R— R mean estimate at a very 
fine timescale (with an arbitrarily small bin size A), which requires no 
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interpolation between the arrival times of two beats. Given the 
proposed parametric model, the nonlinear indices of the HR and 
HRV will be defined as a time-varying function of the parameters 

m = [e(t),goit),gi(0,t),...,gi(P,t),g2iO,0,t),...,g2iQ,Q,t), 
g3(0mt),.:,g3iK,K,K,t)] . 

Parameter Estimation and Model Goodness-of-fit 
measures. A local maximum likelihood method [41] is used 
to estimate the time-varying parameter set ((t). Given a local 
observation interval (t — l,t] of duration /, we consider a subset 
Ufn:n of the R-wavc events, where m= min{A: : Uk>t — l} and 
n = max{A: : Ujc<t}. At each time t, we find the parameter vector 
^(t) that maximizes the local log-likelihood, given the R-wave 
events recorded in the local observation interval: 



A=limsup-log(Lf(OI) 



(18) 



More generally, let us consider an w-dimensional linear system 
in the form yi = Y(t)pi, where Y(t) is a fundamental solution 
matrix with F(0) orthogonal, and {pt} is an orthonormal basis of 
R". Then, the sum of the corresponding n Lyapunov Exponents 
(A/) is minimized, and the orthonormal basis {/?/} is called 
"normal" [54]. One of the key theoretical tools for determining 
LEs is the continuous QR factorization: Y(t)= Q(t) R(t) [55,56] 
where Q(t) is orthogonal and R(t) is upper triangular with positive 
diagonal elements Ru, i=l : n. Therefore we obtain [54—56] : 



L{^{t)\Um:n)= w{t-Uk+x) 



(16) 



log [/■ {uk+ 1 1 , ,Hu^^ 1 m),mt))) ] 



where w(T)=exp(wT) with w = 0.02 , is an exponential 
weighting function for the local likelihood. This value has been 
empirically chosen by considering a range of discrete values 
(w = {0.005,0.01,0.015,0.02,0.025,0.030}), and by choosing the 
optimum according to KS plots goodness-of-fit analysis, as 
described in [41]. We use a Newton-Raphson procedure to 
maximize the local log likelihood in eq. 1 6 and compute the local 
maximum-likelihood estimate of ^{t). Because there is significant 
overlap between adjacent local likelihood intervals, we start the 
Newton-Raphson procedure at t with the previous local 
maximum-likelihood estimate at time / — A in which A define 
how much the local likelihood time interval is shifted to compute 
the next parameter update. We determined the optimal orders 
using the Akaike Information Criterion (AIC) by fitting a 
subset of the data using local likelihood method [41]) as well as the 
Kolmogorov-Smirnov (KS) statistic in the post hoc analysis. More 
in detail, we can compare the AIC scores and choose the 
parameter setup with the minimum AIC value of 
AIC = — 2L + 2 dim (^) where L is the maximized value of the 
likelihood function for the estimated model, and dim(^) is the 
number of parameters in the statistical model. 

It is known from point process theory [41] that the Conditional 
Intensity Function (GIF) P{t) is related to the inter-event 
probability p{t) with a one-to-one relationship: 



m= 



Pit) 



(17) 



The estimated GIF is used to evaluate the goodness-of-fit of the 
proposed heartbeat interval point process probability model, 
which is based on the KS test [41]. 



= lim^log||F(0/^,|| = lim^log||i?(0/^,|| 
= lim \ log||i?,7(0ll, I < i < n 



(19) 



In our case, the matrix Y{t) corresponds to the Jacobian of the 
M-dimensional system of the NARL model parameters, where M 
is the value of the largest model order [39] . Therefore, given the 
NARL model reported in eq. 13 and using eqs. 9 to 12 bringing 
back to the NAR framework, it is possible to obtain an M- 
dimensional state space canonical representation. 

Using the Einstein notation, we obtain: 



liRRiUHtm) = yi(i)ARRi + y2(iJ)ARRiARRj 
+ y3(ij,k) ARRi ARRjARRk 



(20) 



Therefore, the elements of the Jacobian can be computed as 
follows: 



dARRr, 



= yi(i)^p-^ y2(jJ)ARRiSjp + 

y2(iJ)ARRiSip + y^(ij,k)SipARRjARRk (21) 
+ y^(ij\k)SjpARRiARRk 
+ y^(iJ,k)SkpARRiARRj 



and considering the symmetry properties of the NAR kernels, we 
finally obtain: 



dARR„ 



= yi(p) + 2j272(.hP)ARRi 

i 

+ ^J2J2y3(U,P)ARRiARRj 



(22) 



Instantaneous Lyapunov Exponents Estimation 

The Lyapunov Exponent (LE) of a real valued function f(t) 
defined for ^ > 0 is: 



Considering N data samples, we evaluate the Jacobian over the 
time series, and determine the LE by means of the QR 
decomposition: 
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J{n)Q(^n-\) = Q{n)R{n) with ^ = 1 ,2. . .,7V. 



The matrix Y{t) corresponds to thejacobian of this system [39]: 



J{n) = 



0 
0 
0 

0 

dARR, 



1 

0 
0 



dARR2 



0 

1 

0 



dARR} 



Standard and Nonlinear Measures of Heartbeat Dynamics 

In order to perform a comparison analysis with standard and 
nonlinear estimates of heartbeat dynamics, we also calculated the 
standard mean of the RR intervals (Mean RR), the root mean 
square of successive differences of intervals (RMSSD) and the 
number of successive differences of intervals which differ by more 



0 
0 

1 

0 

dARR4 



0 
0 
0 



dARRM 



This decomposition is unique except in the case of zero diagonal 
elements. Then, the LEs 1/ are given by 



tN 



J2 ^^^(^■)" 



(23) 



7 = 0 



where T is the sampling time step. The estimation of the LEs is 
performed at each time / from the corresponding time-varying 
vector of parameters, ^(t). We define the first EE as the 

instantaneous dominant Eyapunov exponent (IDEE). In particu- 
lar, the median IDEE {-Xrr) and its median absolute deviation 
(cT;^^^) were considered as group features. 



than 50 ms (pNN50% expressed as a percentage of the total 
number of heartbeats analyzed) [57]. Referring to morphological 
patterns of HRV, the triangular index is obtained as the integral of 
the histogram (i.e. total number of RR intervals) divided by the 
height of the histogram which depends on the selected bin width 
[57]. Moreover, we performed the estimation of the dominant 
Eyapunov exponent according to the algorithm described by Wolf 
et al. [34] (LdJ and Rosenstein et al. [58] (Lds)- Both algorithms 
are suitably applied to experimental noisy data. Finally, other 
nonlinear measures such as the approximate entropy (ApEn) [59] , 
Sample Entropy [60], and the Detrended Fluctuation Analysis 
(DFA) [61] were evaluated. 
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Figure 1. Simulation results using the Henon equations. (Top) Synthetic interval series from the Henon map, (Middle) the a stepping, and 
(Bottom) IDLE averaged estimate from 100 realizations using NARL model with noise level of /: = 0.001. 
doi:1 0.1 371 /journal. pone.01 05622.g001 
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Experimental Data 

Synthetic Data 

Before reporting the implementation of the synthetic datasets, it 
is important to clarify some important differences from standard 
definitions that are introduced by our methodology. In fact, 
standard nonlinear systems such as the Henon Map and Rossler 
Attractor are intrinsically defined in the continuos time domain, 
whereas our methodology deals with stochastic point processes 
which are a sequence of events. Moreover, additive noise terms 
have to be considered as well. Therefore, starting from the 
canonical equations of each nonlinear system, we slightly modify 
the system equations by adding a noise term. Then, the output of 
the system is taken as an input of an integrate-and-fire system. The 
output of such an integrate-and-fire system constitutes the series 
modeled by the proposed cubic NARL model within the point- 
process framework. 

Henon Map. In order to test the efficiency of the proposed 
cubic NARL model in tracking the complexity of a synthetic 



stochastic series through the IDLE index, we simulated a modified 
version of the well-known chaotic Henon Map as suggested in 
[39]. Such a complex system, in which stochastic terms are also 
considered, is governed by the following differential equations: 

yn+i = bxn + ke{t)] 

(24) 

Xn+\ = 1 - axl + yn + ke{t)] 

The time series yn were taken into account fixing /? = 0.3. The 
term £{t) is an independent identically distributed Gaussian 
random variable with zero mean and standard deviation of 1, 
which is modulated by the coefficient k. The coefficient a is taken 
as a time-varying variable from a=\ to a =1.09 with step size 
0.03. Note that the value a=\.07, in a purely deterministic 
domain, corresponds to the transition between the non-chaotic 
and chaotic regime. A total of 100 realizations the Henon Map 
series were simulated, each of which was comprised of 4000 data 
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Figure 2. Simulation results using the Rossler equations. (Top) Syntlietic interval series from tlie Rossler system, (Middle) the a stepping, and 
(Bottom) IDLE estimates using NARL model. The dotted vertical line indicates the transition between the non-chaotic and chaotic regime. 
doi:10.1371/journal.pone.0105622.g002 
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points with 1000 samples for each of the four a-values. A 
realization of the simulated time series is illustrated in Fig. 1 along 
with the a-values and the corresponding IDLE results. 

Rossler Attractor. As a further validation, we simulated a 
modified version of the well-known chaotic Rossler time series (see 
previous work [19]). Such a complex system, in which stochastic 
terms are also considered, is governed by the following differential 
equations: 



dx 
dt 
d^ 
It 



-z - y + ke(t); 
b -\- z(x-c) + ke(t) 



dy 
dt 



x-\-ay-\-ke{t)] 



(25) 



The time series were implemented with sampling time of 0.01 
using the Runge-Kutta integration and fixing b = 2, and c = 4. The 
term e(t) is an independent identically distributed Gaussian 
random variable with zero mean and standard deviation of 0.01. 
The coefficient a is taken as a time-varying variable from a = 0.35 
to a = 0.45 with step size 0.025. Note that the value a = 0.432, in a 
purely deterministic domain, corresponds to the transition 
between the non-chaotic and chaotic regime. A total of 75000 



data points were generated with 15000 samples for each of the five 
a-values. The simulated time series is illustrated in Fig. 2 along 
with the a-values and the corresponding results on the IDLE. The 
IDLE transitions to positive values reflect the simulated switch to 
chaotic behavior. 

Experimental Data 

In order to validate the proposed algorithms performance as 
related to actual cardiovascular dynamics, we have considered two 
experimental datasets. Since the experimental protocols are fully 
described in [19,41] in this paragraph we provide only briefly 
descriptions of the two datasets. 

Postural Changes. In order to validate the proposed 
algorithms performance as related to actual cardiovascular 
dynamics, we studied the complexity of RR interval series 
recorded from 10 healthy subjects for the study of cardiovascular 
and autonomic regulation undergoing a tilt-table protocol. This 
choice is motivated by the high presence of nonlinearity and non- 
stationarity in such time series. In this study, a subject lying 
horizontally in a supine position is then actively or passively tilted 
to the vertical position. The rest (supine) and upright sessions 
alternated six times with three possible modality of transition: 
active stand-up, slow and fast passive tilt. A single-lead ECG was 




100 200 300 400 500 600 700 

Samples 



Figure 3. Box plots of standard and proposed Lyapunov estimations performed on the Henon map among the <^-stepping with 
noise level of Ar= 0.001. 

doi:1 0.1 371 /journal. pone.01 05622.g003 
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Table 1. IDLE Results from the Henon map Synthetic Dataset. 





a-values 


1 


1.03 


1.06 


1.09 




Noise 0.001 










p-value 




0.1424 + 0.0035 


0.1598 + 0.0041 


0.1814 + 0.0040 


0.1960 + 0.0030 


<10-i2 




0.0004 + 0.0002 


0.0064 + 0.0009 


0.0377 + 0.0019 


0.1772 + 0.0030 


<10-i2 




0.3679 + 0.0195 


0.3530 + 0.0179 


0.4410 + 0.0268 


0.6885 + 0.0349 


<10-i2 




0.0407 + 0.0151 


0.0379 + 0.0112 


0.0636 + 0.0130 


0.0747 + 0.0170 


<10-i2 


Noise ii^O.Ol 










p-value 




2.5488 + 0.0251 


2.5552 + 0.0387 


2.5537 + 0.0280 


2.5573 + 0.0238 


>0.05 




0.3038 + 0.0070 


0.3049 + 0.0077 


0.3037 + 0.0085 


0.3033 + 0.0083 


>0.05 


^RR 


1.2240 + 0.3232 


1.5006 + 0.4380 


1.4822 + 0.5034 


1.4962 + 0.5237 


<5*10^ 


^/^RR 


0.4511 +0.1559 


0.6884 + 0.2008 


0.5780 + 0.2476 


0.7137 + 0.2830 


2<10-5 


Noise ii^O.l 










p-value 




2.6704 + 0.0265 


2.6710 + 0.0283 


2.6756 + 0.0306 


2.6719 + 0.0350 


>0.05 




0.2069 + 0.0997 


0.2161 +0.1221 


0.2320 + 0.1334 


0.2334 + 0.1459 


>0.05 


^RR 


0.1180 + 0.0517 


0.1383 + 0.0565 


0.1319 + 0.0487 


0.1566 + 0.0468 


>0.05 


^^RR 


0.1348 + 0.0376 


0.1398 + 0.0314 


0.1360 + 0.0378 


0.1369 + 0.0247 


>0.05 


Noise /c = 7 










p-value 




0.2949 + 0.0033 


0.2953 + 0.0031 


0.2949 + 0.0031 


0.2947 + 0.0023 


>0.05 




0.0023 + 0.0022 


0.0015 + 0.0015 


0.0011 +0.0011 


0.0017 + 0.0016 


>0.05 


^RR 


-0.1032 + 0.0104 


-0.1018 + 0.0107 


-0.0998 + 0.0098 


-0.0993 + 0.0119 


>0.05 


^?^RR 


0.0318 + 0.0056 


0.0311 +0.0071 


0.0337 + 0.0066 


0.0319 + 0.0071 


>0.05 



doi:1 0.1 371/joumal.pone.01 05622.t001 



KS plots goodness-of-fit analysis, according to [41]. The simulated 
time series along with the resulted IDLE series are shown in Fig. 1 , 
whereas the corresponding box plots are shown in Fig. 3 in terms 
of IDLE median [Xrr) and its median absolute deviation o";^^^. The 
proposed IDLE is able to track the complexity variation at each 
moment in time. As a matter of fact, the IDLE goes increasingly 
high from the non-chaotic behavior to the chaotic one. 
Remarkably, the non-chaos-chaos transition is instantaneously 
detected, although a significant oscillatory dynamics is present in 
the chaotic region. The related IDLE values are reported in 
Table 1 along with standard estimates of the dominant Lyapunov 
exponents. A non-parametric statistical analysis has been per- 
formed in order to quantify the differences between the considered 
a-values for each of the considered noise level. 

Considering the noise level /: = 0.001, the Kruskal-Wallis test 
reveals significant differences (p<10~^^) for both the standard 
Lyapunov estimates and the proposed Xrr and complexity 
variability index cr;^^^. In this case, the Dunn test for multiple 
comparison, which considers a Tukey-Kramer correction, shows 
that each group of standard estimates of Ldj and Ld2 having 
coherent a parameter are different with all the other groups 
(p< 10~^), whereas coherent a-values of Xrr are different with all 
the other groups (p<10~^) except for a = {1 — 1.03} (p>0.05). 

Considering the noise level /: = 0.01, the Kruskal-Wallis test 
reveals significant differences (p<5*10 ^) only for the proposed 
Xrr and complexity variability index O";^^^. In this case, the Dunn 
test for multiple comparison, which considers a Tukey-Kramer 
correction, shows that coherent a=\ values are different with all 
the other groups (p<10-^), with ^ = {1.03,1.06,1.09} equal 
between each other (p>0.05). For noise level k = 0.\ and k=\^ 



continuously recorded for each subject during the study, and the 
RR intervals were extracted using a curve length-based QRS 
detection algorithm [62]. Further details on the experimental 
protocol can be found in [19,41]. The study was conducted at the 
Massachusetts Institute of Technology (MIT) General Clinical 
Research Center (GCRC) and was approved by the MIT 
Institutional Review Board and the GCRC Scientific Advisory 
Committee. Patient records /information was anonymized and de- 
identified prior to analysis. 

Congestive Heart Failure. The second heartbeat dataset 
was constituted from data gathered from Congestive Heart Failure 
(CHF) and reference healthy subjects on a public source: Physionet 
(http://www.physionet.org/) [63]. The RR time series were 
recorded from 14 CHF patients (from BIDMC—CHF Database) 
as well as 16 healthy subjects (from MIT — BIH Normal Sinus 
Rhythm Database). Each RR time series was artifact-free (upon 
human's visual inspection and artifact rejection) and lasted about 
50 min (small segments of the original over longer recordings). 
These recordings have been taken as landmark for studying 
complex heartbeat interval dynamics [3,8]. 

Results 

Instantaneous Complex Dynamics on Synthetic Data 

Henon Map. We performed the IDLE estimation by fitting 
the NARL model on y series from the modified stochastic Henon 
Map time series (see eq. 24). The series were generated a hundred 
times for each of the four considered noise levels 
^ = {0.001,0.01,0.1, 1}. For /:>0.001, a further constrain of 
Xn = 0. le(/) was imposed as x„ < = — 0. 1 or x„ > = — 0.5 in order 
to prevent the system to become unstable. The model orders were 
set as P= 3, 2 = 1, ^ = 1, and a = 0.4 were chosen by preliminary 
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the difference between the a-values are not revealed by standard 
and proposed dominant Lyapunov estimates. 

Rossler System. We performed the IDLE estimation by 
fitting the NARL model on the x series from the modified 
stochastic Rossler time series (see eq. 25). The model orders were 
set as P = 3, Q=\, K=\^ and a = 0.2 were chosen by preliminary 
KS plots goodness-of-fit analysis, according to [41]. The simulated 
time series along with the resulted IDLE series are shown in Fig. 2. 
Clearly, the proposed IDLE is able to track the complexity 
variation at each moment in time. As a matter of fact, the IDLE 
goes increasingly high from the non-chaotic behavior to the 
chaotic one. Remarkably, the non-chaos-chaos transition is 
instantaneously detected, although a significant oscillatory dynam- 
ics is present in the chaotic region. Intervals expressed as median 
± M.A.D. are as follows: -0.0973 + 0.0044 for a = 0.350, 
-0.0813 ± 0.0065 for a = 0.375, -0.0850 ± 0.0077 for 
fl = 0.400, -0.0375 ± 0.0075 for a = 0.425, 0.0033 ± 0.0092 
for a = 0.450. A non-parametric statistical analysis has been 
performed in order to quantify the differences between the 
considered a-values. The Kruskal-Wallis test reveals significant 
differences (p< 10~^) and the Dunn test for multiple comparison, 
which considers a Tukey-Kramer correction, shows that each 
IDLE group having coherent a-values is different with all the other 
group (p<10-^) except for a ={0.375 -0.4} (p>0.05). 

Instantaneous Complex Dynamics on Postural Changes 

Before estimating the IDLE from the experimental datasets, we 
first considered a specific time-domain method [64] for testing the 
presence of nonlinearity in the heartbeat intervals. The null 
hypothesis of the test states that the given time series is linear. In 
the considered recordings, we restricted the test to short-term 
dependence by setting the number of laps M = 8, and a total of 
500 bootstrap replications. Concerning the RR series gathered 
during postural changes, the nonlinearity test shows that the level 
of nonlinearity of the considered RR intervals is statistically 
significant for all the considered subjects but one (see Table 2. As 
also shown in Table 2, the NARL modeling always gives a good 
model fit, with KS distance < 0.0604 in all cases. Specifically, 
concerning the three experimental sessions, i.e. stand-up, slow- tilt, 
and fast-tilt, a decrease of the IDLE with respect to the relative rest 
condition is shown in 25 out of 30 epochs. In particular, in the fast- 
tilt condition the decrease happens for all subjects, and is more 
significant than stand-up and slow-tilt. Group statistics of standard 
and proposed instantaneous measures are shown in Table 3, 
whose inter-subject analysis was performed using a non-paramet- 
ric rank-sum test. Results on the proposed IDLE show a non- 
significant statistical difference between the stand-up epochs and 
their relative rest epochs (/?>0.05) and a significant difference for 
the slow-tilt epochs (/?<0.05). The highest significance was found 
comparing the fast-tilt epochs with their relative rest (/?< 0.001). 
These trends are confirmed by the standard DLL estimation 
according to the Rosenstein et al. [58] technique, whereas the one 
suggested by Wolf et al. [34] did not show such significant 
differences. IDLE dynamics for one representative subject are 
shown in Fig. 4, whereas the averaged IDLEs for all 10 subjects 
are shown in Fig. 5, providing a clear portrayal of how different 
postural stimuli elicit different changes in the dynamic signatures 
of complexity. Concerning other standard and instantaneous 
indices, we report significant differences on three session on /i^^, 
ApEn, and SampEn, whereas RMSSD and pNN50% showed 
significant differences during the slow and fast tilt sessions. 

Using this dataset, we further evaluate the effect of the Laguerre 
parameter a on the IDLE estimates. Tracking for values 
a = {0.1,0. 2,0.3, ...,0.8} from a representative subject undergoing 
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0 200 400 600 800 1000 1200 1400 1600 1800 

Time [s] 

Figure 4. Instantaneous heartbeat statistics computed from a representative subject (subject 1) from the tilt-table protocol. In the 

top panel, the estimated fijiR(t) is superimposed on the recorded RR series. In the bottom panel, the instantaneous averaged IDLE is superimposed 

on the original IDLE. 

doi:1 0.1 371 /journal. pone.01 05622.g004 



postural changes are shown in Fig. 6. Indeed, the IDLE estimates 
are affected by the choice of the Laguerre parameter a. However, 
such a variability is significantly less than the variability of the 
IDLE dynamics within session. As a matter of fact, quantitative 
results reported in Table 4 show that these differences are 
associated to a p-value less than 2e~^ for each experimental 
session. 

Finally, we report further results on the nonlinearity test 
separately performed for each of the experimental session, instead 
of the whole recordings (see Table 2). As a result, under the null 
hypothesis of linearity according to the time-domain method 
described in [64] for testing the presence of nonlinearity in the 
heartbeat intervals, the 49.12% of the resting state (28/57 sessions) 
were associated to a significant p-value <0.05, along with the 
44.4% of the stand-up (8/18 sessions) protocol, the 20% of the 
slow-tilt (4/20 sessions) protocol, the 10.53% of the fast-tilt (2/19 
sessions) protocol. 



Instantaneous Complex Dynamics on CHF patients 

The results of the second experimental dataset (on CHF) are 
shown in Table 5. According to the nonlinearity test, 15 out of 16 
RR time series from the healthy subjects showed significant 
nonlinearity (/? <0.05), whereas in the CHF group, 6 out of 14 RR 
time series failed to reach significance (/?>0.05). The fact that a 
lower degree of nonlinearity was found in the CHF patients 
suggests that pathological conditions might reduce the nonlinearity 
in the heartbeat interval series, which is also consistent with 
previous finding that a healthy heartbeat presents more pro- 
nounced nonlinear dynamics [3,6,8,65]. Table 5 also demon- 
strates that the NARL model well fit both pathological and healthy 
heartbeat series with KS distance < 0.082 in all cases. Results 
averaged among groups are reported in Table 6. We report that 
standard and instantaneous time domain measures are able to 
discern the two groups with high statistical significance 
(p < 5e~^). On the other hand, comparing the standard and 
proposed instantaneous complexity measures, only the DFA-a2 
and the complexity variability a^j^j^ are able to provide significant 
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Figure 5. IDLE dynamics averaged for all 10 subjects. The vertical red lines indicate the transition from the supine to the upright position after 
stand up (top panel), after slow tilt (middle panel), and after fast tilt (bottom panel). 
doi:1 0.1 371 /journal. pone.01 05622.g005 



discrimination capability between the two populations with 
;7<0.05. 

Discussion and Conclusion 

Novelties and Impact of the proposed Methodology 

We presented a novel methodology able to instantaneously 
characterize the complex nonlinear dynamics of a stochastic series 
of events by using the LEs. The proposed approach relies on the 
previous literature for the LEs mathematical definition [38,39] and 
is embedded in a novel IG-based point-process nonlinear 
framework defined through a third-order Wiener- Volterra repre- 
sentation, thus advancing on the previous models [19]. As a 
consequence, the novel instantaneous LEs definition is able to 
provide a reliable complexity measure tool to examine the 
unevenly spaced events at very high temporal resolutions, without 
resorting to any interpolation method. Moreover, goodness of fit 
measures such as KS distance and autocorrelation plots quanti- 
tatively allow to verify the model fit as well as to choose the proper 
model order, which represents another open issue of current 
parametric approaches. 



The effective procedure for the time-varying parameter 
identification is ensured by the combined use of the discrete-time 
Laguerre expansions for the Wiener- Volterra terms and local 
maximum likelihood method. In particular, expanding the 
Volterra terms with the orthonormal Laguerre bases requires a 
reduced number of parameters to retain the information of all the 
past events. The nonlinear regression is further performed on the 
derivative series to better account for nonstationarity [53]. 
Importantly, unlike other methods that might require large sample 
size, our method is potentially useful to perform complexity 
measures in short recordings of the signals of interest. 

Importantly, the proposed measures also allows for the study of 
the complexity variability, i.e., the analysis of complex systems 
referring to the fluctuations in complexity instead of analysis of 
central tendency. Within the proposed framework, it is always 
possible to incorporate physiological covariates (such as respiration 
or blood pressure measures) and produce further instantaneous 
indices from their dynamic cross spectrum and cross bispectrum 
[48]. Unlike other paradigms for estimating nonlinearity indices 
developed in the literature [7,59,66,67], our method is formulated 
within a probabilistic framework specifically developed for point- 
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Table 4. IDLE Variability evaluated through 


a = {0.1,0. 2,0. 3,. ..,0.8} and within each session of the postural changes protocol. 




Session 


Through a-values 


Within Session 


p-Value 


Resting State 


0.0117 + 0.0044 


0.0585 + 0.0213 




Stand-Up 


0.01 25 ± 0.0082 


0.0652 + 0.0178 


2e-' 


Resting State 


0.0170 + 0.0240 


0.0693 + 0.0098 


le-' 


Slow Tilt 


0.0119 + 0.0128 


0.0582 + 0.0113 




Resting State 


0.0167 + 0.0192 


0.0617 + 0.0131 


2e-4 


Fast Tilt 


0.0119 + 0.0128 


0.0500 + 0.0144 




P-values are obtained from the Mann-Whitney test with null hypothesis of equal medians between the two groups. Values are expressed as X = Median(X) + MAD(X). 
doi:1 0.1 371/journal.pone.01 05622.t004 



process observations (e.g. RR intervals), which already produced 
important nonlinear quantifiers for autonomic assessment, based 
on second- and third-order statistics (instantaneous spectrum and 
bispectrum) [19]. Most other nonlinearity indices are derived from 
non-parametric models, whereas our model is purely parametric 
and the analytically derived indices can be evaluated in a dynamic 
and instantaneous fashion. We believe these strengths enable our 
method as a useful tool for assessing nonlinear dynamics of 
heartbeat intervals in a non-stationary environment. 

Study of the Instantaneous Cardiovascular Connplex 
Dynannics 

The novel IDLE index was evaluated in both synthetic and 
experimental heartbeat series. Estimations on the synthetic dataset 
were performed on a stochastic version of the well-known chaotic 
Henon Map and Rossler attractor. The use of such a modified 



version of the Henon Map and Rossler system alongside the 
obtained IDLE results need to be discussed. We are aware that in 
purely deterministic Rossler equations the first EE should be zero 
in the non-chaotic region, whereas it should be increasingly 
positive in the chaotic region. However, the IDEE results show 
slightly negative values in the non-chaotic region. Such a behavior 
may be ascribed to the stochastic input and to the additional 
integrate-and-fire step which affect the estimation of all complexity 
measures, including EEs. To this extent, in order to further 
investigate the effect of noise, results from the Henon Map 
equations were gathered as a function of the noise level. We 
demonstrate that, for small amount of noise (/: = 0.001), standard 
and instantaneous estimates of the dominant Eyapunov exponent 
achieve similar results. However, considering Henon Map 
dynamics with /: = 0.01, the IDEE is exclusively able to discern 
the different behavior of the nonlinear system. Of note, for A: >0.1, 



R 



RTRTRTRTRTR 



a: 
a: 




1500 2000 
Time [s] 



Figure 6. IDLE dynamics averaged for a = {0.1,0.2,0.3,. ..,0.8} from a representative subject undergoing postural changes. Resting state 
(R) sessions during supine position alternate with upriglit session after tilt (T). The bold line and gray area indicate the IDLE median and MAD, 
respectively. Averaged values for all subjects are shown in Table 4. 
doi:1 0.1 371 /journal.pone.01 05622.g006 
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Table 5. Results from the CHF-Healthy Experimental Dataset. 



Subject 


Group 




p-value 


KS dist. 






01 


CHF 


995.4 + 26.4 


>0.05 


0.0445 


0.2268 


0.1123 


03 


CHF 


910.25 + 28.9 


>0.05 


0.0552 


0.2165 


0.1257 


04 


CHF 


603.09 + 22.7 


<0.02 


0.0456 


0.0676 


0.0896 


05 


CHF 


655.6+13.3 


>0.05 


0.0297 


-0.01 14 


0.0507 


06 


CHF 


637.4+15.9 


<0.001 


0.0802 


0.0757 


0.0659 


07 


CHF 


778.0 + 7.3 


>0.05 


0.0363 


0.0762 


0.0793 


08 


CHF 


800.1 + 14.1 


<0.01 


0.0357 


-0.0965 


0.0327 


09 


CHF 


602.5 + 5.4 


>0.05 


0.0305 


— 0.0121 


0.0555 


10 


CHF 


486.9 + 6.9 


<0.001 


0.0329 


-0.0622 


0.0444 


1 1 


CHF 


685.1 +16.0 


<0.02 


0.0354 


-0.0014 


0.0578 


12 


CHF 


722.8 + 27.2 


<0.001 


0.0326 


-0.0798 


0.0380 


13 


CHF 


619.7 + 5.1 


<0.001 


0.0386 


0.0550 


0.0655 


14 


CHF 


837.7 + 23.4 


>0.05 


0.0367 


0.0041 


0.0613 


15 


CHF 


652.15 + 20.6 


< 0.005 


0.0265 


-0.0322 


0.0550 


16265 


healthy 


1023.9 + 38.9 


<0.001 


0.0527 


0.0162 


0.0438 


16272 


healthy 


924.7 + 30.6 


<0.001 


0.0539 


-0.0304 


0.0557 


16273 


healthy 


1 046.2 + 68.7 


<0.001 


0.0764 


0.0998 


0.0543 


16420 


healthy 


849.9 + 39.2 


<0.001 


0.0394 


0.0107 


0.0410 


16483 


healthy 


818.5 + 24.9 


<0.001 


0.0336 


-0.0254 


0.0396 


16539 


healthy 


831.5 + 47.9 


<0.001 


0.0592 


0.0736 


0.0485 


16773 


healthy 


1238.9 + 74.7 


<0.001 


0.0819 


0.0842 


0.0613 


16786 


healthy 


945.1 +37.0 


<0.001 


0.0503 


0.0226 


0.0527 


16795 


healthy 


889.4 + 61.8 


<0.001 


0.0442 


-0.0214 


0.0462 


17052 


healthy 


939.8 + 33.3 


>0.05 


0.0487 


0.0682 


0.0819 


17453 


healthy 


816.6 + 31.9 


<0.001 


0.0416 


0.0310 


0.0466 


18177 


healthy 


639.3 + 25.7 


<0.001 


0.0267 


0.0058 


0.0608 


18184 


healthy 


831.2 + 31.7 


<0.001 


0.0392 


-0.0373 


0.0438 


19090 


healthy 


993.54 + 41.0 


<0.001 


0.0565 


-0.0101 


0.0393 


19140 


healthy 


849.1 +49.9 


<0.001 


0.0353 


0.0459 


0.0513 


19830 


healthy 


821.2 + 29.4 


<0.001 


0.0452 


-0.0795 


0.0332 



P-values are obtained from the nonlinearity test. 
doi:1 0.1 371/journal.pone.01 05622.t005 



the noise has an amplitude comparable with the output of the 
system, thus destroying the different behaviors among the a-values. 

The use of the Laguerre expansion of the Wiener- Volterra 
kernels was also investigated through experimental analysis. As the 
zero-order Laguerre basis is an exponential function, the IDLE 
estimates present a mild dependence on the a value of the 
Laguerre functions. Nevertheless, we demonstrated that the actual 
information needed to characterize the experimental sessions, i.e., 
the variability within each session, is significantly higher than the 
variability among all the a values. Anyway, using the hereby 
proposed approach we clearly demonstrate the ability of the IDLE 
in tracking the system complexity in an instantaneous fashion. The 
IDLE, in fact, becomes higher when the simulated system switches 
from non-chaotic to chaotic behavior (see Fig. 2). In all 
applications, an IG probability model was used as a stochastic 
version of the widely-applied deterministic integrate-and-fire 
models used to simulate heartbeats. Regarding the experimental 
datasets, we demonstrated that our approach is useful in 
characterizing the inherent nonlinearity of the cardiovascular 
system. For the first time, tracking complexity by instantaneous 



Lyapunov Exponents was performed and evaluated during 
postural changes. During the resting condition the cardiovascular 
and autonomic nervous system are more sensitive to the initial 
conditions (positive IDLE), whereas a more regular dynamics 
(negative IDLE values) appear during the tilt phases (see Fig. 4). 
These results are in agreement with previous findings that complex 
vagally-driven dynamics are blunted under sympathetic drive [68] 
and with more recent reports on loss of complexity during states of 
arousal [20]). Our instantaneous measures also confirm that the 
instantaneous complexity reflects instantaneous autonomic ner- 
vous system (ANS) control on the cardiovascular dynamics. We 
have shown that tracking ANS complexity on healthy subjects 
undergoing postural changes not only confirms previous results 
[69,70], but further improves sympathovagal assessment as elicited 
by different dynamic gravitational stimuli. 

Our experimental findings on the nonlinearity test performed 
on each experimental session (resting state, stand-up, slow tilt, and 
fast tilt) suggest that loss of instantaneous heartbeat complexity as a 
function of velocity of the postural changes is reasonably due to 
changes in the nonlinearity of the cardiovascular system (instead 
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Table 6. Group Statistics of Standard and Instantaneous Heartbeat Dynamics Measures from CHF-Healthy Experimental Dataset. 







CHF (n = 14) 


Healthy (n = 16) 


p-value 


Standard and Instantaneous Time Domain Measures of HRV 


Mean RR (ms) 


669.73 ± 68.66 


855.74 + 56.14 


< 0.002 


RMSSD 


0.0121 +0.0036 


0.0432 + 0.0145 


<4£?-4 


pNN50% 


0.2357 + 0.2246 


21.5406+15.4908 


<1e-4 


HRV Triangular Index 


2.9551 +0.5769 


2.5628 + 0.3593 


>0.05 




671.55 + 69.6 


864.7 + 53.3 


<4e-3 




8.31 +2.2 


24.7 + 7.0 


<5e-4 


Standard and Instantaneous Nonlinear Measures of HRV 


ApEn 


1.2130 + 0.1032 


1.2177 + 0.1066 


>0.05 


SampEn 


1.5670 + 0.2690 


1.4092 + 0.1522 


>0.05 


DFA-ai 


0.8498 + 0.2191 


1.0820 + 0.1467 


>0.05 


DFA-a2 


1.1552 + 0.1335 


0.9286 + 0.0544 


<0.05 


Ld, 


0.0167 + 0.0025 


0.0165 + 0.0012 


>0.05 


Ld, 


0.0029 + 0.0008 


0.0033 + 0.0005 


>0.05 


^RR 


0.0014 + 0.0649 


0.0135 + 0.0368 


>0.05 




0.0595 + 0.0120 


0.0476 + 0.0066 


<0.05 


P-values are obtained from the Mann-Whitney test with null hypothesis of equal medians between the CHF and healthy subject groups. Values are expressed as 



X = Median(X) + MAD(X). 
doi:1 0.1 371/journal.pone.01 05622.t006 



of, for example, changes on the noise properties). Of note, we have 
previously reported that the standard HRV indices defined in the 
time and frequency domain are unable to distinguish the three 
possible modality of transition through different p-values [19]. 

Moreover, the novel complexity features {^rr; c;.^^} give 
important information in the complexity evaluation, also useful in 
distinguishing heartbeat dynamics coming from patients with CHF 
and healthy subjects. We found that pathological heartbeat 
dynamics are associated with increased complexity variability, 
providing an unique measure of complexity able to discern the 
CHF and healthy populations. Of note, some of the standard 
measures defined in the time domain have similar p-value than the 
point-process measures, as well as DFA-a2 shows similar 
performances than the IDLE median absolute deviation. Never- 
theless, we point out that we aimed at showing the performances of 
novel instantaneous measures of complexity based on Lyapunov 
exponents, providing novel insights on the complexity character- 
ization of stochastic time-varying discrete point-process systems. In 
other words, although other HRV-based measures are able to 
discern CHF from healthy subjects, no other measures have been 
proposed to characterize the time-varying complexity behavior 
occurring in pathological vs. a healthy cardiovascular system. In 
particular, while confirming the results reported in the current 
literature (i.e., several measures of complexity are not able to 
characterize the CHF and healthy subjects groups), we show a 
novel key complexity behavior through the proposed complexity 
variability framework. These findings can be linked to the current 
literature whereby cardiovascular disorders affect complexity and 
variability, and may lead to serious pathological events such as 
heart failure [71]. 

Concerning other preliminary applications, the proposed IDLE 
methodology has been revealed as a powerful tool also in tracking 
the instantaneous complexity during loss of consciousness induced 



by anesthetic drugs [46]. Looking at the overall results shown 
using actual heartbeat dynamics data, it seems that the median 
IDLE is sensitive to changes in ANS regulation induced by 
orthostatic stress, while IDLE median absolute deviation is 
sensitive to changes in ANS regulation induced by CHF. 
However, this conclusion cannot be made without speculation at 
this time as the data used for the comparison between CHF 
patients and healthy subjects is related to long-term ECG 
monitoring during unstructured activity, whereas the data from 
the tilt-table experimental dataset is structured. Therefore, the 
observed sensitivity of the IDLE measures could be due to 
different physiological behavior occurring in CHF subjects or to 
differences related to the kind of experimental protocol. Future 
works are related to pursue this direction in further investigating 
the potential of these high-order nonlinear models in producing 
new real-time measures for the underlying complexity of 
physiological systems, and to the investigation of the instantaneous 
complexity, along with the baroreflex sensitivity and respiratory 
sinus arrhythmia, during postural changes in CHF subjects, thus 
allowing some conclusions on the complex physiological behavior 
of the cardiovascular system in CHF subjects. 
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